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Abstract 

We obtain a modified version of the Onsager regression relation for the expectation values of 
quantum-mechanical operators in pure quantum states of isolated many-body quantum systems. 
We use the insights gained from this relation to show that high-temperature time correlation 
functions in many-body quantum systems can be controllably computed without complete diago- 
nalization of the Hamiltonians, using instead the direct integration of the Schroedinger equation for 
randomly sampled pure states. This method is also applicable to quantum quenches and other sit- 
uations describable by time-dependent many-body Hamiltonians. The method implies exponential 
reduction of the computer memory requirement in comparison with the complete diagonalization. 
We illustrate the method by numerically computing infinite-temperature correlation functions for 
translationally invariant Heisenberg chains of up to 29 spins 1/2. Thereby, we also test the spin 
diffusion hypothesis and find it in a satisfactory agreement with the numerical results. Both the 
derivation of the modified regression relation and the justification of the computational method 
are based on the notion of quantum typicality. 
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In 1931, Onsager came up with the profound insight that "the average regression of fluctu- 
ations will obey the same laws as the corresponding macroscopic irreversible process"-. This 
statement known as the "Onsager regression hypothesis" (ORH) became the cornerstone of 
the linear response theory. From today's perspective, the ORH is equivalent-^ to the high- 
temperature limit of the fluctuation-dissipation theorem^. In this paper, we show that a 
modified version of ORH holds for the expectation values of quantum-mechanical operators, 
when a many-body system is in a pure quantum state. We also present an efficient method 
for computing high-temperature linear response characteristics of many-particle quantum 
systems using the time evolution of a single pure state. 

There exists a class of nonperturbative problems, such as nuclear spin-spin relaxation 
in solids-, where the relaxation or correlation functions in translationally-invariant systems 
need to be computed at high temperatures. Despite the progress in the approximate meth- 
ods, e.g.- 1 ^, and numerical techniques^ - — , the above kind of problems generally resist con- 
trollable solutions, leaving the complete diagonalization of quantum Hamiltonians as the 
only way to obtain controllable results. The sizes of the systems treatable by complete di- 
agonalization are severely limited by the computer memory requirement that scales as N 2 , 
where N is the number of quantum states in the system. The memory requirement for the 
controllable-accuracy algorithm proposed in this work scales at most a iV(logiV) 2 . 

In recent years, it was realized that, given the exponentially large number N of quan- 
tum states in a many-particle system, many observable properties of such a system can be 
obtained by sampling one suitably chosen pure quantum state, or a wave function — the 
so-called "quantum typicality"— - —. In particular, Refs.— - — applied the notion of quantum 
typicality to the relaxation and fluctuation phenomena, but on the numerical side these 
investigations dealt so far only with the systems that were sufficiently small to perform the 
complete diagonalization of the Hamiltonians. 

In this paper, we report a conceptual and a computational results, which are both con- 
nected to the notion of typicality but, otherwise, only indirectly connected to each other. 
The conceptual result is that the expectation values of quantum-mechanical operators in a 
pure quantum state obey the usual regression relation but with the amplitude of fluctuations 
exponentially reduced in comparison with the classical case (see Eq. ([8|) below). The compu- 
tational result is that the high-temperature time correlation functions can be controllably 
computed on the basis of Eq.flH]) without complete diagonalization of the Hamiltonians, us- 
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ing instead the direct integration of the Schroedinger equation for randomly sampled pure 
states. As an example, we obtain infinite-temperature correlation functions for translation- 
ally invariant Heisenberg chains of up to 29 spins 1/2, thereby also testing the spin diffusion 
hypothesys. To the best of our knowledge, none of the complete diagonalization studies of 
the Heisenberg spin-1/2 chains conducted so far has reached the above size. We note here 
that pure quantum states were used in Refs.— ^ in the context of other numerical methods. 

Below, in order to be specific, we consider a lattice of N s interacting spins 1/2 with the 
total number of quantum states N = 2 s ^> 1 and the Hamiltonian H. We adopt the 
following conventions: (i) Analytical formulas are presented only in the leading order in 
1/N. (ii) Wave functions without time argument and operators with time argument imply 
the Heisenberg representation of quantum mechanics. The opposite implies the Schrodinger 
representation, (iii) h — 1. 

We now focus on some observable quantity, e.g. total spin polarization, characterized by 
operator A = A mn , which has zero average value at the infinite-temperature equilibrium, i.e. 
Tr{/i} = 0. The Onsager regression relation for this quantity near the infinite temperature 
equilibrium has the following form: 

Tr{i(t)p ncq } = ^Tr{i(t)i(0)}. (1) 

where A(t) = e Mt Ae~ int , and p neq = -^exp(aA) with a being a small constant. The 
right-hand side (RHS) of Eq.flT]) represents the equilibrium time correlation function of A, 
while the left-hand side (LHS) is the relaxation function of A corresponding to the initial 
nonequilibrium density matrix p neq . 

Quantum typicality investigation of Ref.~ implied that 



(^ neq |i(t)|^ neq ) = Tr \A(t)p 

ncq 



o 



1 



(2) 



where |\I/neq) is a wave function that "samples" p neq . 

Now we obtain a complementary relation on the fluctuation side. It involves the wave 
function \^ e q) representing the infinite temperature equilibrium and defined as a random 
vector in the Hilbert space of the system. |^/ e q) can be generated in any orthonormal basis 
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{|0fc)} as follows: 

N 

l*eq) = (3) 
fc=l 

where a k are the quantum amplitudes, whose absolute values are selected from the proba- 
bility distribution^^ 

P(\a k \ 2 ) = Nexp(-N\a k \ 2 ) (4) 

and the phases are chosen randomly in the interval [0,2-71"). In the following, we use bar 
above an expression to indicate the Hilbert-space average over all possible choices on |^/ e q)- 
Now we consider the correlation function for the time series of the expectation value 
{^ eq \A(t')\^ eq ) in the time interval [— T, T + t): 

C(t,T) = ±j\>(y cq \A(t + t>)\y eq )(y eq \A(t')\y eq ). (5) 

IdM- we derive the following relation: 

C(t,T) = ±Tr{A(t)A(0)}+A 1 , (6) 

where 



A? « 2 ^ TN , f *b ( [ Tr {ife)i(O)}] 2 + Tr [A(t - t 2 )A(0)) Tr [A(t + f a )i(0)}) . 

~ TV ~ 2 (7) 
For large enough T, the correction Ai in Eq.(j6]) is much smaller than the principal term as 
long as Tr |a(£)A(0) j *~^°°> 0. In particular, if Tr |A(t)A(0)| decays at large t faster than 
|£|-o.5 then, f or i ar g e enough T, the integral in Eq.^7]) becomes independent of T, and, as 
a result, Ai = 0(y^r/T)Ti |^ 2 | /N 2 , where r is the characteristic timescale for the decay 
of the expression under the integral. If Tr |vi(t)A(0) j decays asymptotically as \t\~ u with 
< v < 0.5, then, according to Eq.(j7J), A x still remains small, but its prefactor scales as 
0(T- U ). 

Eq.Q together with Eqs. (TTH2T) implies the modified version of ORH for the expectation 
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value of the operator A in a pure quantum state: 



lim {y neq \A(t)\V neq ) = a lim NC{t,T), (8) 

N— >oo ' T— >oc,N— >oo 

where, in the RHS, the limit iV — > oo is taken first, which in practical terms means that T 
should be much smaller than the inverse spacing of the energy levels in the system as the 
above limits are taken^. 

From the viewpoint of practical computing, the implications of the above find- 
ings are two- fold: (i) As already implicit in Eq.(|2]), and explicit in Eqs.( )6|8|) . a 
single realization of {9 neq \A(t) \^ n eq) is exponentially more accurate in approximating 
Tr |y4(t)y4(0)| than the corresponding single classical relaxation process in approximat- 
ing classical correlation function 23 . (\I/ neq |yi(t)|\I/ neq ) decays into the equilibrium statisti- 
cal noise (^ eq |A(t)|^ cq ), which, according to Eq.® has root-mean-squared (rms) ampli- 



tude a/C(0, T) m yTr{A 2 }/N, which is by factor of y/N smaller than the rms amplitude 



Tr{A 2 }/N expected for the classical noise or the noise of continuously monitored macro- 
scopic quantum observable at infinite temperature^ 3 - 1 ^. This noise suppression is due to the 
fact that the time evolution of a single pure state contains the superposition of A" independent 
"noises" associated with each of the basis states^. The statistical noise of neq \A(t)\^ ncq ) 
can be suppressed further by averaging over many pure-state evolutions, (ii) In principle, as 
we show below, the direct evaluation of C(t,T) can also be used to obtain Tr jyt(t)A(0) j, 
but this procedure does not take advantage of the above-mentioned quantum parallelism 
and hence is less efficient. 

Although the evaluation of (\I/ neq |v4(t)|\l/ neq ) is a very efficient method to obtain 
Tr |v4(t)v4(0) j, an even more efficient method is to use typicality to sample this trace 
directly on the basis of the following relation derived in^: 

(*^\A(t)A(0)\**d = {Mt)A(0)} + A 2 , (9) 

where 

~%=±T:{A(t)A(0)A(t)A(0)}. (10) 
That the second term in the RHS of Eq. (Q is much smaller than the first one can be shown 
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FIG. 1: Intermediate dynamic structure factor I n (t) of the Heisenberg chain of 20 spins 1/2. 
The calculation based on the exact trace formula (jlip is compared with the approximations given 
by Eqs. (|2)6I9|) as indicated in the legend. The initial agrement between the black and the green 
lines also demonstrates the validity of the regression relation ([8]). All calculations are based on 
the complete diagonalization of the Hamiltonian T~L. Each of the three approximate calculations 
is done with a single pure state. In the case of Eq.([2]), a = 0.083 corresponding to the initial 
polarization equal to approximately 4 percent of the maximum polarization. In the case of Eq.©, 
T = 4200/ J. As expected theoretically the approximation based on Eq.© gives the most accurate 
agreement with the exact result. (Note the logarithmic vertical scale.) The accuracy of all three 
approximations can be improved by averaging over more pure states. 



,/Trji 4 ) 

by estimating their ratio at t = as >ry r \ 2 i ~ "7KF- ^ ne statistical accuracy of computing 



Tr{A*} 



Tr ^A(t)A(0)j with the help of Eq.(j9]) is thus better by factor 1/a in comparison with 
Eq.©. 

In FigJUwe demonstrate the relationships (EJ EJ |9]) by computing the intermediate dynamic 
structure factor I % (t) for the Heisenberg chain of 20 spins 1/2 using complete diagonalization. 
Thereby we also demonstrate the regression relation (JS}. The Hamiltonian of this chain is 
H = JJ2i Si • Si + i with periodic boundary condition. Here J is the coupling constant, and 
Si is the spin operator on the ith chain site. Such a chain admits periodic spin modulations 
with wave numbers q = 27in/N s , where n is an integer number taking values 0, 1, N s — 1. 
For a given wave number q, the intermediate dynamic structure factor is denned as 



I q (t) = Tr \A {q} (t) A {q} (0) 



(11) 



where A [q} = T /m cos ( ( l m ) S m- 

Now, we proceed with showing that, for the spin systems too large to be treated by 
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complete diagonalization, it is still possible to controllably compute infinite temperature 
correlation functions by evaluating the LHS of Eq.(j9]) with the help of the direct integration 
of the Schrodinger equation. 

We compute the time evolution of pure states on the basis of the time-discretized version 
of the Schrodinger equation. We use the fourth-order Runge-Kutta routine based on the 
following equation: 

\V(t + At)) = \V(t)) + h) + \v 2 ) + \v 3 ) + \vi), (12) 

where At is the discretization time step, and \vi), \v 2 ), \v 3 ), |t> 4 ) are unnormalized Hilbert- 
space vectors computed as follows: \vi) = — %H\^{t)) At, \v2) = — ^iH\v\)At, \vz) = 
— ~H-L\v2)At, and ^4) = — \iTV\V3) At. Given the linearity of the Schrodinger equation, 
Eq.( ll2p is equivalent to the simple 4th order power-series expansion of the time-evolution 
operator at each discretization time step. We used At = 0.01/ J 

The above routine requires only storing in the memory the vectors \^f) and \vi) and the 
non-zero elements of the Hamiltonian H. Although the Hamiltonian is an N x iV matrix, it is 
very sparse for many-particle systems with only two-particle interactions when represented 
in a "local" basis, where each basis function is factorizable in terms of the wave functions 
of individual particles. For N s spins 1/2, a possible local basis is the one where the z- 
projections of all spins are quantized. In this basis, the number of the nonzero entries of the 
Hamiltonian matrix is of the order of N x iV 2 for the systems with long-range interactions, 
or iV x N s for the short range interactions. Thus the overall memory required for the direct 
propagation of the Schroedinger equation scales at most as iV(log N) 2 , i.e. it is exponentially 
smaller than the memory required for the complete diagonalization, which scales as iV 2 . One 
can take advantage of this memory reduction only when the operator of interest, A, is also 
sparse in the local basis, but this is normally the case in physical contexts. In fact, in many 
cases, including the calculations of I q (t), it is possible simply to use the eigenbasis of A as 
the local basis. 

We verify the accuracy of the direct integration method in two ways. For small spin 
clusters, we compare the wave functions obtained by propagating the same initial state using 
either complete diagonalization or the direct integration method. As shown in Fig. |2fa), the 
overlap between the two wave functions remains extremely close to 1 over the time interval 
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FIG. 2: Tests of numerical accuracy of the direct integration of the Schrodinger equation, (a) 
Overlap between two initially identical wave functions for the Heisenberg chain of 20 spins 1/2, 
| x I / Ex(i)) and \^RKi(t)), computed using, respectively, the complete diagonalization and the direct 
integration, (b) Overlap between two initially identical wave functions for the Heisenberg chain 
of 29 spins 1/2, and \^2(t)), both computed using the direct integration method with two 

respective discretization time steps At\ = 0.01/ J and At 2 = 0.02/ J. The noisiness of the line 
originates from the accumulated machine rounding errors. 

required to compute I n (t) for 20 spins 1/2 in Fig. [TJ For larger systems, we compare two 
wave functions and \^2(t)) obtained by propagating the same initial wave function 

using the direct integration method with two different discretization time steps Ati and At 2 
such that At2 = 2At\. We then verify that that their overlap {^i(t)\^2{t)) is sufficiently 
close to 1. An example of such a test for 29-spin Heisenberg chain is shown in Fig. 12(b). 

We note here that the same direct integration algorithm can be used to compute the 
imaginary-time evolution associated with the expression exp(— "H/3/2)|\l/ eq ), where /3 is the 
inverse temperature, thereby generating equilibrium wave function corresponding to tem- 
perature 1//3. This wave function can then be used to compute the linear response charac- 
teristics at temperature 1/(3. 

Our numerical procedure for computing (\l/ eq |A(t)A(0)|\l/ eq ) is based on propagating 
two wave functions using the direct integration method. One of them is \^ eq (t)) = 
exp(— i7it)\^/ cq (0)) , where |\I/ eq (0)) is given by Eq.([3]). The other one is |$(t)) = 
exp(-i%t)|$(0)), where |$(0)) = A\^ cq (0)) (i.e. |$(0)) is unnormalized) . The quantity 
of interest (^ cq |y4(t)A(0)|^ eq ) is then evaluated as (^ cq (t)\A\$(t)). 

Now we exemplify the direct integration method by computing the intermediate dynamic 
structure factors I q (t) with q = 2tc/N s for Heisenberg chains of sizes up to N s = 29. By 
doing this calculation, we also test the spin diffusion hypothesis, which stipulates that, for 
sufficiently small values of q, Ig(t) ~ e~ Dq2t , where D is the diffusion coefficient. 

Our results presented as plots (a) in Figj3] indicate that, in every case, I q {t) shows the 
initial tendency to decay exponentially, but then the behavior universally starts exhibiting 
oscillations. Plots (b) in the same figure further indicate that the nearly exponential parts of 
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FIG. 3: Intermediate dynamic structure factors I q (t) for Heisenberg chains of N s spins 1/2 
computed as the LHS of Eq.Q by propagating a single pure state with the help of the direct 
integration of the Schrodinger equation. The values of N s are indicated in the plot legend. In 
each case, q = 2ir/N s . The horizontal axis is: (a) f{Jt) = Jt, (b) f(Jt) = £i q 2 t, and (c) 
f(Jt) = £2 q 2 (1 + 0.1 Inq) tint, where £1 and £2 are arbitrary scaling parameters. The vertical axes 
for (b) and (c) are displaced for better visibility. Plots (a) represent the original calculation results. 
Plots(b) test the scaling expected for spin diffusion. Plots (c) test the empirical scaling reported 
for classical spins in Ref. 2 -. The numerical accuracy test for the pure state evolution in the case 
of the 29-spin chain is given in Fig. [2](b). 

I q {t) exhibit satisfactory g 2 -scaling, while plots 2(c) show that the scaling g 2 (l + 0.1 lng)lnt 
reported in the numerical studies of classical spins 2 ^ works even better. 

To summarize, we obtained the modified Onsager regression relation ([8]) for a pure quan- 
tum state. We also find that the direct computation of the LHS of Eq.(j9]) is the most 
efficient way to obtain equilibrium time correlation functions with controllable accuracy. 
We have directly tested only the high-temperature limit but the method itself can also be 
used at finite temperatures. We further note that the direct integration of the Schrodinger 
equation in combination with the random sampling of pure states can also be used for the 
efficient computing of quantum quenches and other situations describable by time-dependent 
many-body Hamiltonians. 
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knowledges the hospitality of the Kavli Institute for Theoretical Physics at the University 
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National Science foundation under Grant No. NSF PHY11-25915. The numerical part of 
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SUPPLEMENTAL MATERIAL 

Note: Equation numbers, figure numbers and citation numbers appearing in this supple- 
mental material start with letter "S". Equation, figure and citation numbers without "S" 
refer to the text of the main article. 



I. DERIVATION OF EQS.(6,7) 

To obtain Eq.(6), the correlation function (5) can be represented in the energy eigenbasis 
as follows: 

i r T 

C(t,T) = —l dt' a\ala v a m A qv A nm ^-^^-^-^ (SI) 

— m,n,p,q 



Due to the random phases of complex amplitudes a n , the Hilbert-space average a* q a* n a p a m is 



not zero only when the four amplitudes form conjugate pairs such as a*a q a^a n . The main 
contribution to the RHS of Eq.f lSlj) comes from the 2-pair terms where the indices q and n 



are different from each other, in which case, according to Eq.(4), a*a q a^a n = (|a g | 2 ) 2 = 1/N 2 . 
The number of the remaining terms, where the two pair indices coincide, is smaller by factor 
1/N, and therefore, their contribution to the RHS of Eq. (IS II) is suppressed by the same factor 
as long as the set of the eigenvalues of A does not have a small subset of the anomalously 
large members. Limiting ourselves to the terms with different pair indices, we obtain 



a*a* n a p a m = 1/N (S gp 5 nm + 5 qm 5 np ) (S2) 

The term 5 qp 8 nm in the parentheses does not make contribution to the average because of the 
condition Tr{A} = 0. The summation over the term S qm 5 np with the subsequent integration 
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over t' gives 

W7T) = ^- 2 J2e t{n)t A mn A nm = ^Tr{i(t)i(0)} , (S3) 

m,n 

which is the principal term in the RHS of Eq.(6). 



The deviation of C(t,T) obtained for a particular realization of ^f eq from C(t,T) can be 

The explicit form of this variance 



estimated from the variance A 2 = C(t,T) — C(t,T) 
obtained using Eq.( tS2|) is 

A 2 — — !— f T rlt' C T rlf" A A A,, A p*(e g -ep+efc-ei)t p i(e q -e p +e n -em)t' p i(e k -ei+e u -e v )t" 
1 — 4T 2 }— r p u,h }_'Y Ulh J ^-qp- ri -nm- ri -kl- ri -uv e e c 

,m,n,k,l,u,v _iyi (Pqp^nrn ^qm^np)^kl^uv 

5 kv 5 ul )} . (S4) 



The evaluation of a*a^a* k a^apa m aia v relies on the same considerations as those that led 
to Eq. (lS2p . Namely, the nonzero contribution to the sum in the RHS of Eq. (lS4p comes 
from the terms, where the eight amplitudes are organized into conjugate pairs such as 



a* q a q a* n a n a* k akO* u a u . The main contribution to the RHS of Eq. (lS4j) comes from the 4-pair 
terms which have all four indices q, n, k and / different from each other. In this case, according 



to Eq.(4), a*a q a* n a n a\a]~a* u a u = (|a g | 2 ) 4 = 1/iV 4 . Now we observe that the number of all 



possible combinations of four different conjugate pairs in expression a q ^* rl a\a lf u a p a m aia v is 
relatively large. However, we find by inspection that most of these combinations eventually 
give zero contribution to the RHS of Eq. (lS4p because of the condition Tr{A} = 0. Below, 
we only include those combinations that give nonzero contributions: 

a v — Jjz[ &qrn&np&kv&ul + &ql&nv$kp$um + $qv&nl$km&up (S5) 
OqmVnl3kv3up ^qv^np^km^ul ^qm^nv^kp^ul 

Oql5 n p5k v 5 um ~\- fiqlfinv3km3up ^qv^nl^kp^um 

+ {omitted terms} ]. 



After substituting Eq. (IS5j) into Eq. (IS4p most of the resulting terms contain integral of 
the type: 

1 r T r T 

' dt'dt" e iJlt ' e"""", (S6) 



AT 2 j_ r j_ r 

where rj is some combination of energies e n . The integration region for the above integral is 



12 



shown in Fig JSlf a). Since our primary goal is to find the order of magnitude of the leading 
terms in the limit of large T, we adopt the following relatively crude approximation, which 
significantly simplifies the analytic expressions. Namely, we rotate the integration region as 
shown in Fig JSlf b) and simultaneously change the integration variables to i' x = (£' +t")/y/2 
and t x = (t" — t') j y/2. Finally, we substitute t 2 = £ x y/2, thereby replacing expression ( IS6l) 
with 

dt 2 e~ ivt2 . (S7) 



2TV2 J^tV2 

Among the nine terms in the square brackets in Eq. flS50 , not all terms contribute equally 
to the RHS of Eq.( lS4l) . The main contribution comes from the first, the second and the third 
terms. The contribution from the first term is then canceled exactly by the contribution of 
the term containing 1/N 4 in the square brackets in Eq. (lS4p . What remains after that is: 

Af « 2 ^ Tm dt 2 ( [Tr {i(t 2 )i(0)}] 2 + Tr [A(t - i 2 )i(0)} Tr [A(t + t 2 )i(0)}) . 

~ TV " (S8) 
The two terms in the above integral originate, in their respective order, from the second and 
the third terms in Eq. (1S5j) . Each of the remaining six terms in Eq. (lS5j) would contribute 
to the RHS of Eq.flSi} a term of the following type: Tr ^A(t)A{0)A(t + t 2 )A(t 2 )}. For 
t 2 ~ T ^> t, the operators separated by t 2 should become uncorrelated, which would imply 
that the above term becomes equal to Tr |a(£)A(0)| , which is smaller than the terms 
included in the estimate flS8|) by factor 1/N. 

In general, the two terms in Eq.f ]S8j) are of the same order of magnitude, when t is of the 
order of the characteristic decay time of the correlation function Tr |A(i)A(0) |. However, 
if Tr{i(t)i(0)} decays asymptotically faster than |t| 0,5 , then the second term in Eq. (IS8j) 
becomes much smaller than the first one as t increases. 

Since the above derivation is done in the leading order in 1/N, it implicitly assumes that 
T is much smaller than the inverse level spacing of the system. If one is to consider the 
limit T — > oo at fixed N, then the contributions from the corrections of higher order in 1/N 
would need to be examined. 
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FIG. SI: (Color online) Figures illustrating the rotation of coordinates and the rotation of the 
integration region, (a) Original coordinates (t',t") and the original integration region, (b) Rotated 
coordinates^,^) and rotated integration region. 

II. DERIVATION OF EQS.(9,10) 

The derivation of Eq.(9) is based on the following general relationships: If V is a Hermi- 
tian operator, then we can use the eigenbasis of V to define | \l/ eq ) in Eq. (3) and then also use 
the probability distribution (4) to obtain (^ eq |V r |^ r e q) = J2 m \ a m\ 2 V mrn = jjTr^Vj. The 
variance with respect to the above average is: A| = J2 m (l a m| 2 — l/N) 2 V^ m = ^Tr |^ 2 |- 
(This formula is different from the closely related Eq.(2) of Ref. [17], because, as noted in 
[22], we allow statistical noise for the normalization of |^ e q)-) 111 the context of Eq.(9), we 
can choose V = A(t)A(0) and thereby obtain Eq.(lO). 



III. COMPARISON WITH THE RELAXATION OF CLASSICAL SPINS 

In order to illustrate the computational advantage of the quantum parallelism associated 
with the direct evaluation of (\l/ neq |A(t)|\l/ neq ), let us consider what it takes to calculate 
the high-temperature relaxation function of the z-component of the total magnetization of 
classical spin systems numerically. 

Let us assume that the system investigated numerically is a chain of N s classical spins 
governed by the Hamiltonian 

— ^ JijSixSjx + JfjSiySjy + JijSi Z Sj Z , (S9) 

i<j 

where (S ix , S iy , S iz ) are the projections of classical spin vectors having length 
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yj Sf x + Sf y + Sf z = 1, and J?-, J?-, J?- are the anisotropic coupling constants. 

We are interested in the linear response regime. Hence the initial polarization of the 
spin system should be much smaller than N s . Let us take it to be equal to 0.1N S . At 
the same time, the rms magnetization noise level for this system at the infinite temperature 
equilibrium is yj N s (S^ 2 ) = y/N s /3. Thus, in a single numerical realization of the relaxation 
process, one can only obtain a reasonable accuracy in the range between 0.1 N s and roughly 
2^/N s /3. For 1000 classical spins, this means that a single run of a relaxation process will 
only be statistically accurate in the range between 100 and 40. The standard way to improve 
the quality of the calculation is to average over n independent realizations of the relaxation 
process. The statistical noise in this case decreases very slowly — only as 1/ y/no. 

If one computes (^/ neq |A(t)|\l/ neq ) for a single realization of |\&neq) for N s spins 1/2 with the 
initial 10 percent polarization, i.e. (^ f n eq|^.(0)|^ f n eq) = 0.1 x ^N s , then this single run would 
only be affected by the statistical noise associated with (^/ eq |74(t)|^/ eq ) when (\l/ neq |A(t)|\I/ ncq ) 

decays to the values of the order of 2 ^1/™ • ^ or ^ s = 1000, sucn a noise would be too 
small to be of any practical concern. However, for a system of 20 spins, the equilibrium 
noise does have observable effect on {^ neq \A(t)\^/ ncq }, as can be seen in Fig. 1. As in the 
classical case, this noise can be further suppressed by factor l/y/no after averaging over n 
independent realizations of (^ neq |A(t) |\l/ neq ). 



IV. MONITORING A 

The fluctuations of (\l/ eq |A(t)|\l/ e q) should be distinguished from the measured equilibrium 
noise, when a macroscopic quantum observable A is continuously monitored [S1-S4]. Let us 
again assume that A represents the total z-component of the magnetization of the spin 
system. In this case, the rms amplitude of the monitored infinite-temperature spin noise 
is expected to be \J ±Tr{A 2 } = y / 'n s (S^ 2 ) [24], while, as implied by Eq.(6), the rms 
amplitude of the quantum mechanical expectation value is suppressed by the extra factor 
1 / \/N . As mentioned in the main article, this suppression is the consequence of the fact 
that the noise of the quantum-mechanical expectation value is, in a sense, averaged over 
the quantum superposition of iV independent "noises" associated with the time evolution of 
each of the basis states in the Schroedinger representation. 

Despite the exponentially large difference in the amplitudes, the monitored noise and 
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the noise of {^ cq \A(t)\^ cq ) should have identical normalized time correlation functions in 
macroscopic systems. If the total magnetization is the only variable being continuously 
monitored in the system (or one of a small number of variables), then the disturbance of 
the dynamics of an individual spin in the bulk of the sample should be minimal. Since the 
fluctuations of the total magnetization come mostly from spins in the bulk, this means that 
the correlation function of the monitored noise should still be proportional to the quantum 
mechanical trace Tr |A(t)A(0)| computed for a perfectly isolated quantum system. 
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